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Abstract 


Two Reynolds-averaged Navier-Stokes computer codes - one unstructured and one structured - are 
applied to two workshop cases (for the 3rd Workshop on CFD Uncertainty Analysis, held at Instituto 
Superior Tecnico, Lisbon, in October 2008) for the purpose of uncertainty analysis. The Spalart- 
Allmaras turbulence model is employed. The first case uses the method of manufactured solution and 
is intended as a verification case. In other words, the CFD solution is expected to approach the exact 
solution as the grid is refined. The second case is a validation case (comparison against experiment), 
for which modeling errors inherent in the turbulence model and errors/uncertainty in the experiment 
may prevent close agreement. The results from the two computer codes are also compared. This 
exercise verifies that the codes are consistent both with the exact manufactured solution and with 
each other. In terms of order property, both codes behave as expected for the manufactured solution. 
For the backward facing step, CFD uncertainty on the finest grid is computed and is generally very 
low for both codes (whose results are nearly identical). Agreement with experiment is good at some 
locations for particular variables, but there are also many areas where the CFD and experimental 
uncertainties do not overlap. 


1 Description of the Computer Codes 

FUN3D [ 1—3] is a finite-volume RANS solver (either compressible or incompressible equations can 
be solved) in which the flow variables are generally stored at the vertices of the mesh. FUN3D 
solves the equations on mixed element grids, including tetrahedra, pyramids, prisms, and hexahedra 
and also has a two-dimensional path. It employs an implicit upwind algorithm in which the inviscid 
fluxes are obtained with a flux-splitting scheme. At interfaces delimiting neighboring control vol- 
umes, the inviscid fluxes are computed using an approximate Riemann solver based on the values 
on either side of the interface. For second-order accuracy, interface values are obtained by extrapo- 
lation using gradients computed at the mesh vertices using an unweighted least-squares technique. 
Limiting of the reconstructed values may be employed for flows with strong shocks. For all re- 
sults presented in this paper, the convective flux scheme used is Roe’s flux difference splitting [4], 
For tetrahedral meshes, the full viscous fluxes are discretized using a finite-volume formulation in 
which the required velocity gradients on the dual faces are computed using the Green-Gauss theo- 
rem. On tetrahedral meshes this is equivalent to a Galerkin type approximation. For non-tetrahedral 
meshes, edge-based gradients are combined with Green-Gauss gradients, which improves the h- 
ellipticity of the operator, and the complete viscous stresses are evaluated. The solution at each 
time-step is updated with a backward Euler time-differencing scheme. At each time step, the linear 
system of equations is approximately solved with either a multi-color point-implicit procedure or an 
implicit-line relaxation scheme [5], Local time-step scaling is employed to accelerate convergence 
to steady-state. FUN3D is able to solve the RANS flow equations, either coupled or uncoupled 
with the Spalart-Allmaras [6] (SA) one-equation turbulence model. The Menter SST Model [7] is 
also available for uncoupled solutions. In this paper all computations are uncoupled and use the SA 
model. By default, the turbulence advection terms are discretized using first-order upwinding. 

An emerging capability in FUN3D is the ability to solve with a cell-centered discretization in- 
stead of node-centered. Although this capability is not fully production-ready at this time, some 
preliminary results will be shown for the manufactured solution. Also, in these cell-centered so- 
lutions, the turbulence advection terms are discretized with a conservative second-order accurate 
scheme. 

CFL3D is a structured-grid upwind multi-zone CFD code that solves the generalized thin-layer 
or full Navier-Stokes equations. [8] For all results in this paper, the full Navier-Stokes equations 
have been employed. The code can use point-matched, patched, or overset grids and employs local 
time-step scaling, grid sequencing and multigrid to accelerate convergence to steady stage. A time- 
accurate mode is available, and the code can employ low-Mach number preconditioning for accuracy 
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in computing low-speed steady-state flows. CFL3D is a cell-centered finite-volume method. It uses 
third-order upwind-biased spatial differencing on the convective and pressure terms, and second- 
order differencing on the viscous terms. Roe’s flux difference-splitting (FDS) method [4] is used 
to obtain fluxes at the cell faces. The solution is advanced in time with an implicit approximate 
factorization method. A wide variety of eddy-viscosity turbulence models are available in the code, 
including nonlinear models. Only SA has been used in this study. Cross-derivative terms are ignored 
in the turbulence model. The turbulence models are solved uncoupled from the mean flow equations, 
and, unless otherwise specified by the user, are solved with first-order upwind turbulence advection 
terms. 


2 Discretization Uncertainty Estimation 


The current procedure for the estimation of discretization uncertainty is based on the grid conver- 
gence index (GCI) from Celik et al. [9], The GCI on the fine grid is given by: 
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and (j > i is the quantity of interest on the finest grid, </> 2 is the quantity of interest on the next-finest 
grid, T 21 is the ratio of cell spacing from one grid to the next, and p is the order of accuracy of the 
method. 

In the current paper, the grids are all structured (solved as mixed-element type hexahedra in the 
unstructured code), and successively coarser levels are constructed by removing every other grid 
point in each coordinate direction. Thus, r2i = 2 for all results herein. When three grid levels are 
available and ?’2i = r.32 = 2, the order p is computed from: 
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where £32 = fa — fa and £ 2 i = fa ~~ fa- When the quantity £ 32/ £21 < 0- the convergence is said to 
be “oscillatory,” and when l^l < \£ 2 i \ the sequence is “divergent.” In either case, eq. (3) no longer 
applies. 

The discretization uncertainty, U, of the solution on the fine grid is defined by: 


U = GCI 21 \fa\ (4) 

However, as Eca and Hoekstra [10] point out, if p < 1, uncertainty estimates tend to be over- 
conservative, and, for p much higher than the theoretical order of the method, the uncertainty esti- 
mates can be unreliable. Therefore, we adopt the following methodology for estimating the uncer- 
tainty (some of these ideas are from Eca and Hoekstra): 

• For 0.95 < p < 3.05, U = GCI 21 \fa | 

• For 0 < p < 0.95: U = min(GCI 21 \fa\ , 1.25A M ) 

• Forp > 3.05: U = max(GCI* 21 \fa\ , 1.25A m) 

where Am is the maximum difference in absolute value between the fa on any of the three grid 
levels used to determine p, and GCI* 21 is the same as eq. (1) except that p is taken to be 3: 
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For oscillatory convergence, the uncertainty is taken as: U = GCI** 21 \<pi\, where a presumed 
order of convergence p = 2 is assumed, and a higher factor of safety is employed (2.5 as opposed to 
1.25): 
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For divergence (|£ 3 2 | < |e 2 i|), the uncertainty is taken as U = 3A m- 


( 6 ) 


3 Manufactured Solution 

The manufactured solution used is the same as that described as MSI in Eca et al. [11] for flow on 
a square domain 0.5L < x < L and 0 < y < 0.5 L, namely: 

u/u ref = erf(rj) (7) 

= <8) 

— = 0.5(n(2a: - x 2 + 0.25)ln(4y 3 - 3 y 2 + 1.25) (9) 
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where rj = a y/x, cr = 4, ry, = o v y/x, a v = 2.5cr, and v max = 10 3 ^. The Reynolds number 

Re = u re fLlv = 10 6 . The variable v is the turbulence variable from the SA model, which is 

related to the eddy viscosity v t via v t = vf v i, where f v i = (v / v) 3 / [(y / v) 3 + 357.911]. 

Note that FUN3D can be run using incompressible equations, but CFL3D is only a compressible 
code. When solving the compressible equations, an additional exact solution variable of constant 
total temperature everywhere in the domain was specified, based on an assumed freestream Mach 
number of 0.2. Also for the compressible equations it was necessary to scale u, v, and p appropri- 
ately. 

Although not shown, the codes were also run using the MS2 and MS4 exact solutions, but 
in these cases the turbulence variable could be driven negative very near the wall over a portion 
of the domain. This behavior is possibly a result of the fact that these turbulent exact solutions 
have an asymptotic behavior near the wall that is dramatically different from the generally accepted 
behavior [12, 13] of v t oc y 3 . In fact, for MS2: v t oc y 8 , and for MS4: v t oc y 16 . For MSI 
the behavior is more reasonable: ex y 1 . In Eca et al. [11] special handling of the turbulence 

manufactured source term helped to avoid numerical difficulties; in FUN3D and CFL3D no such 
special handling was performed. 

It should also be pointed out here that the SA model used in FUN3D and CFL3D has minor 
differences from that employed by Eca et al. In particular, Eca et al. leave out the term f t 2 , which 
is a part of the model as specified in Spalart and Allmaras [6]: 
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( 12 ) 


ft 2 = c t3 exp (-c t4 x 2 ) 

Leaving out the ft> term has a minor influence on the solution behavior. We include the f t 2 term 
into our source terms defining the SA exact solution. Also, in FUN3D and CFL3D, an important 
computational limit is placed on the variable r (which goes into the computation of /,,,). The variable 
r should remain positive to keep the solution well-behaved: 
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and So is a very small positive number. It was noticed for MSI that if S is not clipped for eq. (13), 
the forcing term becomes unbounded at the locations where r passes through zero. 

For the manufactured solution, the exact solution was specified along the left, top, and right 
boundaries for FUN3D (node-centered) and in two rows of ghost-cell centers for CFL3D. For 
FUN3D (cell-centered), the exact solution and gradient for the viscous terms were specified along 
these boundaries. The bottom wall used standard no-slip boundary conditions, with dp/dn = 0 
(and, for compressible flow, either dT/dn = 0 in CFL3D or T = constant in FUN3D). The turbu- 
lence variable v at the wall was set to zero. 

Two sets of stretched grids were used. The first set, termed “sgl,” had a fine grid size of 577 x 
577. Spacing was uniform in the x-direction (Ax = 0.00086806), and was stretched in the y- 
direction. The first grid point off the bottom wall was at y = 4.3403 x 10 -5 . A total of 5 nested 
grid levels were employed: 577 x 577, 289 x 289, 145 x 145, 73 x 73, and 37 x 37. The second 
grid set, termed “sg2,” had a fine grid size of 289 x 1153. Spacing was uniform in the x-direction 
(Ax = 0.00173611), and was stretched in the y-direction. The first grid point off the bottom wall 
was at y = 1 x 10 -6 . A total of 5 nested grid levels were employed: 289 x 1153, 145 x 577, 
73 x 289, 37 x 145, and 19 x 73. 

To determine the discretization error, the exact manufactured source terms were added to the 
Navier-Stokes equations, and the codes were converged to near machine-zero (10 -14 ) residuals. 
Various integral and point quantities were compared with the exact solution. 

Fig. 1 shows error in drag (from the exact solution c ( j = 16ln(2)p / ^/tt = 6.25706270620 x 
10~ 6 ) for results on sgl using the two codes with both first- and second-order turbulence advection 
terms. The “(cc)” indicates that the cell-centered version of FUN3D was used for one of the runs. 
In this plot, N represents the total number of degrees of freedom (e.g., number of cell-centers or 
grid points), so \/l/N is a measure of the average grid spacing h. All results for this integral 
quantity converge to the exact solution with second-order behavior on the finest grid levels. Note 
that when using the compressible equations, it is necessary to remove the effect of variable density 
when post-processing drag to compare with the exact solution specified for incompressible flow. 

It is important to note that although the absolute value of errors in drag coefficient are monoton- 
ically decreasing, the actual drag levels do not exhibit monotonic convergence until the finest three 
grid levels. This can be seen in Fig. 2. The Cd, GCI, and uncertainty on the finest grid level are listed 
in Table 1 . 

Figs. 3 and 4 show LI -norm errors for various quantities, using second-order and first-order 
discretization for turbulence advection, respectively. In Fig. 3, convergence is second-order for all 
quantities, as expected. For Fig. 4, convergence degrades to first-order. Results on the sg2 grid (with 
finer normal spacing) are given in Fig. 5. Again, when the turbulence advection terms are discretized 
second-order, the global error norms for all quantities converge second-order. It should be noted in 
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Table 1. The c,j and its uncertainty on the finest grid for MSI on sgl 


Code 

Turbulence advection order 

Cd 

approx GCI, % 

U(c d ) 

CFF3D 

2nd 

0.644531 x 10 -5 

3.05 

0.196577 x KT b 

CFF3D 

1st 

0.638512 x 10" 5 

4.26 

0.271690 x 10" 6 

FUN3D 

2nd 

0.644844 x 10“ 5 

5.09 

0.327975 x KT 6 

FUN3D 

1st 

0.637872 x 10" 5 

7.67 

0.489037 x 10 -6 


practice, however, that for most aerodynamic problems of interest, little practical influence of first- 
order treatment of the turbulence advection terms has been found using typical grids. 

Fig. 6 shows the convergence of various flow field quantities at a specific location in the domain 
as a function of h both grid sets. Both codes are approaching the exact solution, with the sg2 series 
of grids generally converging in a “smoother” fashion than the sgl series. Note that in some cases, 
the convergence is oscillatory (non-monotonic), even on the finest grid levels. 


4 Backward Facing Step Validation Case 

The backward facing step case was run on the same series of hexahedral grids in both codes. Only 
the node-centered formulation of FUN3D was employed, and both codes were run using first-order 
discretization for turbulence advection. The finest grid consisted of 255, 266 grid points (253, 952 
cells), with minimum spacing at the bottom and top walls of 1.5 x 10 -4 . This spacing yielded an 
approximate average y + at the first cell-center of about 0.25 on the finest grid. The grid was not 
clustered with viscous spacing along the back face of the step, however: minimum spacing there 
was 0.0208333. Successively coarser grids were constructed by removing every -other point in each 
coordinate direction. Fig. 7 shows the grid 3-levels-coarser than the finest grid (with only 3968 
cells). The grid clustering from the upstream near-wall region continued into the shear layer, with 
some spreading. 

All solid walls were solved as no-slip walls. In CFL3D these were treated as adiabatic, whereas 
in FUN3D they were given a constant temperature equal to the freestream adiabatic temperature. At 
the upstream boundary, the velocity and turbulence were specified according to the manufactured 
solution given by the workshop organizers. At the downstream boundary, pressure was specified 
( p/Pref = 1-00149), and all other variables were extrapolated from the interior of the domain. 

Convergence of various surface integral quantities with decreasing grid spacing h (= y/l/N) 
is shown in Fig. 8. These drag coefficients are defined in the standard way, as a force divided by 
loo S re f, where q ^ = l/2p ao u' 2 00 . In previous uncertainty workshops, e.g., Eca et al. [14], the drag 
forces were nondimensionalized without the 1/2 factor in the denominator. As a result, the current 
drag results are a factor of 2 larger than earlier reported results. Both codes generally go to similar 
results as the grid is refined, although for c,l v on the bottom wall it is not clear whether the results 
will cross or stay together if the grid could be refined further. 

Estimates of the uncertainty are also shown for the finest two grid levels. These estimates are 
computed using the methodology described earlier, using the grid level in question plus 2 coarser 
levels. For Cd, v on the bottom wall, CFL3D is well-behaved with an apparent (approximately) 
second-order convergence, so its uncertainty levels are small; FUN3D on the other hand is con- 
verging at less than first-order according to the finest three grids, so its uncertainty for this quantity 
is higher. For Cd, v on the top wall, CFF3D exhibits oscillatory convergence, whereas FUN3D is 
monotonically converging. Finally, both codes show well-behaved convergence and hence reason- 
able uncertainty levels for Cd,p . In this case, FUN3D appears closer to the grid-converged result on 
any given grid, so its uncertainty levels are smaller than those of CFF3D. 
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Table 2. Various quantities and their uncertainty on the finest grid for backward facing step 


Code 

quantity 

fine grid value 

approx GCI, % 

U{c d ) 

CFL3D 

Cd tV bottom wall 

0.0523858 

0.07 

0.342654 x 10" 4 


CcL,v t0 P wall 

0.0944691 

0.57 

0.538355 x KT 3 


c p step 

0.2010955 

2.45 

0.491621 x 10" 2 


reattachment 

6.01137 

0.26 

0.0154735 

FUN3D 

Cd, v bottom wall 

0.0523199 

2.69 

0.140690 x 10" 2 


Cd,v top wall 

0.0940830 

0.65 

0.609768 x KT 3 


c p step 

0.2061919 

0.52 

0.106451 x 10" 2 


reattachment 

6.06012 

0.70 

0.0422834 


The reattachment location is shown in Fig. 9. It appears that both codes will yield a reattachment 
point on an infinitely-refined grid of approximately 6.02. 

Table 2 gives values for the various integrated quantities as well as the reattachment point on the 
fine grid, along with the corresponding computed GCI and uncertainty. 

Figs. 10, 11 and 12 show comparisons of wall pressure coefficient and skin friction coefficient 
with experiment. Fine grid results along with error bars for both CFD and experiment are shown. Re- 
sults for the two codes are almost indistinguishable, but agreement with experiment is only marginal 
along the bottom wall. In particular, C p is too low for CFD for x < 3, and Cf is too low for 
CFD both over the first half of the bubble as well as downstream of reattachment. Interestingly, the 
location for reattachment (where C / passes through zero) appears to agree very well with the exper- 
imental location. Although not everywhere completely within the experimental error tolerance, the 
computed top wall pressure coefficients agree very well with experiment. 

Finally, Figs. 13-21 show comparisons of profiles at 3 locations downstream of the step with 
experiment. Fine grid results along with error bars for both CFD and experiment are shown. Again, 
results for the two codes are almost identical. Generally speaking, computed results compare very 
well with experiment, with the notable exception of n-velocity at x - 6//, which is underpredicted 
significantly by the CFD. 


5 Conclusions 

Two RANS computer codes - FUN3D and CFL3D - have been analyzed using the manufactured 
solution MSI, then applied to a backward facing step computation and compared with experiment. 
For the manufactured solution on the finest grids, the expected asymptotic behavior was observed: 
discretization error was either first- or second-order accurate, depending on the treatment of the tur- 
bulence advection term. Both codes converged appropriately to the exact solution, although some- 
times with oscillatory-convergence behavior. For the backward facing step, both codes were again 
consistent with each other, in the sense that the solutions generally approached nearly identical re- 
sults as the grid was refined. Comparisons with experiment included uncertainty estimates for both 
CFD and experiment. 
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Figure 2. Drag coefficients for MSI on sgl as a function of h, showing oscillatory convergence on 
the coarsest two grid levels for three of the methods. 
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Figure 3. Ll-norm of errors in u, v, p, and v for MSI on sgl using FUN3D (cell-centered) and 
CFL3D, both with second-order discretization for turbulence advection. 
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Figure 4. Ll-norm of errors in u, v, p, and v for MSI on sgl using FUN3D (node-centered) and 
CFL3D, both with first-order discretization for turbulence advection. 
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Figure 5. Ll-norm of errors in u, v, p, and v for MSI on sg2 using FUN3D (cell-centered) and 
CFL3D, both with second-order discretization for turbulence advection. 
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Figure 6. Convergence of various flow field quantities at x = 0.6, y = 0.001 for MSI on sgl and 
sg2 using FUN3D (cell-centered) and CFL3D, both with second-order discretization for turbulence 
advection. 
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Figure 7. Hexahedral grid for the backward facing step case (showing only every 8th point in each 
coordinate direction for clarity). 
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Figure 8. Convergence of various surface integral quantities, including uncertainty estimates on the 
finest 2 grid levels. 
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Figure 9. Convergence of reattachment location behind step, including uncertainty estimates on the 
finest 2 grid levels. 
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Figure 10. Surface pressure coefficient along the bottom wall, including uncertainty error bars for 
both CFD and experiment. 
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Figure 11. Surface skin friction coefficient along the bottom wall, including uncertainty error bars 
for both CFD and experiment. 
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Figure 12. Surface pressure coefficient along the top wall, including uncertainty error bars for both 
CFD and experiment. 
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Figure 13. Profiles of u/u re f at x = 1 H, including uncertainty error bars for both CFD and 
experiment. 
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Figure 14. Profiles of v/u re f at x = 1 H, including uncertainty error bars for both CFD and experi- 
ment. 
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Figure 15. Profiles of u'v' /v%. e f at x - 1 //, including uncertainty error bars for both CFD and 
experiment. 
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Figure 16. Profiles of u/u re f at x = 4 H, including uncertainty error bars for both CFD and 
experiment. 
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Figure 17. Profiles of v/u re f at x = 4 H, including uncertainty error bars for both CFD and experi- 
ment. 
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Figure 18. Profiles of u'v' /v%. e f at x = 4 H, including uncertainty error bars for both CFD and 
experiment. 
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Figure 19. Profiles of u/u re f at x = 6 H, including uncertainty error bars for both CFD and 
experiment. 
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Figure 20. Profiles of v/u re f at x = 6 // , including uncertainty error bars for both CFD and experi- 
ment. 


27 



Figure 21. Profiles of u'v' /u%. e f at x = 6 // , including uncertainty error bars for both CFD and 
experiment. 
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